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Abstract: A non-perturbative determination of the axial current improvement coeffi¬ 
cient ca is performed with two flavors of dynamical improved Wilson fermions and plaque- 
tte gauge action. The improvement condition is formulated with Schrodinger functional 
boundary conditions and enforced at constant physical volume. Large sensitivity is ob¬ 
tained by using two different pseudo-scalar states in the PCAC relation. We estimate the 
resulting correction to Fps at (3 = 5.2 to be around 10%. 
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1. Introduction 

The approach of lattice observables to their continuum limit can be understood in terms 
of Symanzik’s effective low energy theory [1, 2], When applied to QCD with Wilson 
fermions [3], close to the continuum limit, it predicts that scaling violations are domi¬ 
nated by terms linear in the lattice spacing a. These effects can be removed by adding 
a single term to the action, with properly tuned coefficient Cg W [4], Linear effects in the 
lattice spacing can also be systematically removed in (on-shell) matrix elements of com¬ 
posite operators [5,6]. In particular, 0(a)-improvement of the axial current requires to add 
one dimension four operator, with coefficient ca (see eq. (|2.3|), below). In the quenched 
approximation, the coefficients c sw and ca have been determined non-perturbatively in the 
relevant range of bare coupling (or lattice spacings) in [7]. The improvement conditions, 
which determine the improvement coefficients, were derived from the chiral symmetry of 
the continuum limit. More precisely, the PC AC relation was required to hold at finite 
lattice spacing [5-7]. As will be detailed in sect. 2.2, the PCAC-relation can be consid¬ 
ered with different external states. In [5-7], finite volume states were chosen, formulated 
in the framework of the Schrodinger functional. Later, ca was also estimated from the 
PCAC relation in large volume [8,9] at a couple of values of the lattice spacing. Around 
a ~ 0.1 fm, the results for ca obtained from the finite volume definition [7] differ quite 
significantly from those obtained in large volume in [8,9]. At smaller lattice spacing, the 
difference decreases. 

For the interpretation of this difference, one should keep in mind that beyond pertur¬ 
bation theory the improvement coefficients themselves are affected by 0(a) ambiguities. 
In some detail this has been discussed and demonstrated numerically in [10]. The 0(a) 
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ambiguity simply corresponds to the fact that the improved theory is treated up to 0(a 2 ) 
effects. While this forbids a unique definition of the improved theory, the O(o) ambiguities 
can be made to disappear smoothly if the improvement condition is evaluated with all 
physical scales kept fixed, e.g. in units of r o ~ 0.5 fm [11], while only the lattice spacing is 
varied [10]. We call this the constant physics condition. At the same time one has to take 
care that the improvement conditions are imposed using low energy states with E <C a -1 
since 0(a)-improvement is valid for those only. So far, the methods of [8,9] have not yet 
been implemented such as to satisfy these conditions. 

In [12], two improvement conditions for ca were studied, which are easily generalized to 
respect the above criteria. They are formulated in finite volume with Schrodinger functional 
boundary conditions, which furthermore helps to render the numerical evaluation feasible 
in full QCD. We will discuss the improvement conditions briefly in sect. 2.2 and choose 
one of them to compute ca in the Nf = 2 theory, where c sw is known from [13,14], The 
knowledge of ca is crucial in order to be able to determine the pseudoscalar decay constants, 
but also in order to compute renormalized quark masses starting from the PCAC masses, 
as has been done e.g. in [15,16]. 


2. Strategy and techniques 

Before going into the details of our strategy and techniques let us comment again on the 
constant physics condition. In the Schrodinger functional and neglecting the choice of 
the quark mass for a moment, the relevant point is the following. We need to know how 
the lattice spacing depends on /3 in order to determine the latter such that a certain L/a 
corresponds to a prescribed value of L/r o- This has to be enforced only with a rather 
moderate precision, since (sticking with ro as the reference scale) a relative error A in the 
estimate of r^/a translates into an error of the improvement constant which is proportional 
to a /Lx A. Thus even if A varies a bit in the considered range of lattice spacings, this is 
quite irrelevant, in particular if A is a smooth function of the lattice spacing. 

In the remainder of this section we will state in more detail how the constant physics 
condition is implemented. We will also discuss the methods in [7,12] to determine ca, with 
emphasis on the one we finally used. 


2.1 Constant physics condition 

With two degenerate flavors of light quarks, the theory has two bare parameters (3 and m®. 
The bare quark mass mo controls the physical quark mass and the bare coupling determines 
the lattice spacing, defined at vanishing quark mass (for a more precise discussion, which 
however is of little relevance in this context, see [5,17]). Non-perturbative estimates of 


tr 0 (P) 


[r 0 /a\ (5.2) 

[r-o/a] (/3) 


( 2 . 1 ) 


are available in a limited range of f3 [18,19]. In [20], the results of [18,19] were extrapolated 
to zero quark mass. Taking directly these values for [ro/a](f3), we have the points with 
error bars in Fig. [|. From those we roughly estimated the location of the filled points, using 
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Figure 1: The evolution of the lattice spacing a with the inverse bare gauge coupling (3 from 
perturbation theory in the lattice scheme and large volume data [18] for the scale ro- The filled 
points correspond to our “scaled” simulations. 


the perturbative dependence of the lattice spacing on (3 as a guideline. For our action this 
is given to three loops in [21], which builds on various steps carried out in [22-28]. Applied 
as a pure expansion in the bare coupling (no tadpole improvement), one has 


a(gp) 

«((^) 2 ) 


= e [*o 2 (so) 2 ^ 2bo [gl/{%) 2 ] blpb ° 
q = 0.4529(1), g 0 < g' 0 . 


[ 1 + Q bo — (So)"] + O ((ff'o) 4 ) ] > 


( 2 . 2 ) 


The evolution of the lattice spacing relative to our reference point at (^q) 2 = 6/5.2 is then 
expressed by the function t(/3) = a(6//3)/a(6/5.2), which is plotted as a thick line in the 
graph. It confirms that the filled points are very reasonable choices. Note that other forms 
of applying bare perturbation theory (differing from eq. ( |2.2| ) in the go-term) would give 
somewhat different results, but since we are interested in a rather limited range in g$ this 
does not matter much. Later we will show that systematic uncertainties in ca introduced 
by this approximate scale setting are negligible. 

Finally, we keep the PCAC mass approximately constant. A precise definition will be 
given in the following. 

2.2 Improvement conditions for the axial current 

In this section we discuss criteria for the choice of the improvement condition. With the 
isovector axial current A“(x) = ^{x)\T ar )^^ip{x) and the corresponding pseudo-scalar 
density P a (x) = 'tp(x)^T a 'y 5 'i/j(x) we define the improved axial current [5] 

(^i)S(z) = ^(*) + + d *^) pa ( x ) > ( 2 - 3 ) 
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where d^ (<9*) denotes the forward (backward) lattice derivative. In the following we will 
consider quark masses derived from the PC AC relation 


m(x\ a, /3) 


{alljdf, + d*)(Ai)“(x)\f3) 
2(a\P a (x)\/3) 


(2.4) 


Since this mass is obtained from an operator identity, it is independent of the states |a) 
and |/3) as well as the insertion point x up to cutoff-effects. Enforcing this independence at 
finite lattice spacing leads to possible definitions of improvement conditions [12]. Inserting 
the expression for the improved current (^]|) in the previous equation, the quark mass can 
be written as m = r + cicas + 0 (a 2 ) with 


and 


r(x;a,/3) 


■s(x; a, /3 ) 


(ulhidy + d*)A“(x)\fi) 
2(a\P a (x)\(3) 
(a\d^d;p a (x)\p) 
2(a\P a (x)\/3) ' 


(2.5) 

( 2 . 6 ) 


If we now consider two sets of external states and two insertion points, the improvement 
condition m{x\ a, /3) = m(y, 7 , 5) yields 


^ _ Ar _ 1 r(x;a,/3)-r{y\'y,8) ^ ^ 

aAs a s(x] a, /3) — s(y; 7 , 6) 

and therefore the sensitivity to ca is given by aAs. 

Once a reasonably large sensitivity is achieved, all improvement conditions at constant 
physics are equally valid in the sense that 0 (a) effects are removed in on-shell quantities. 
However, the way in which higher-order lattice artifacts are modified will depend on the 
concrete choice of the improvement condition. In particular, if states with energy not so 
far from the cutoff are involved, large 0 (a 2 ) effects might be introduced. 

We now specialize to Schrodinger functional boundary conditions, introduced in [29,30] 
and recall the definition of the relevant correlation functions. One considers QCD in a 
finite volume L 3 xl with Dirichlet boundary conditions in time and periodic boundary 
conditions in space. More precisely, the fermionic fields are periodic up to a phase 6. 
By taking functional derivatives with respect to fermionic boundary source fields one can 
define correlation functions involving the quark fields £ at xo = 0 and at xq = T. In 
this work we use 

= ( 2 . 8 ) 

X 

/p(zo;w) = -^f]T(p“a)0“M> (2.9) 

X 

and f 1 (u\u) = -^(0 ,a {J)0 a {u)) ( 2 . 10 ) 

with the pseudo-scalar operator 

O a (u) = a 6 ^ C(x) 7 5T a ^(x - y)C(y) (2.11) 

x >y 
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at the xq = 0 boundary and the corresponding operator 0' a (u>') at the upper boundary 
of the SF cylinder. These operators depend on spatial trial “wave functions” u> and u/, 
respectively. 

The Schrodinger Functional version of eqs. (|2.5|) and f[2.(i[) is then given by 


and 


r(x o;w) 
s(x 0 ;w) 


5 (dp + <9 q)/a(R);^) 
2/p(x 0 ;w) 
dodo/ppEoi^) 
2/p(x 0 ;w) 


( 2 . 12 ) 

(2.13) 


To determine ca in the JVf = 0 theory [7], Ar and As were originally defined through 
a variation of the periodicity angle 9 of the fermion fields, while keeping xq = T/2 and 
uj = const fixed. For this method the sensitivity aAs is quite low when L > 0.8fm, 
T = 2L. In addition, with dynamical fermions different values of 9 would require separate 
simulations. We therefore consider this method as too expensive and disregard it in the 
following. In the quenched approximation two alternatives have been explored in [12]. 

Requiring the quark mass to be independent of xo (for fixed 9 and u> = const ) is 
technically easy to implement. However, also in this case the sensitivity is small unless large 
values of 9 are used. Moreover, the contribution of excited states is not well controlled, 
because one insertion point must be rather close to a boundary to achieve a sufficiently 
large sensitivity. Thus energies which are not far removed from the cutoff may contribute. 

Secondly, variations of the wave function u have been considered. Ideally, one would 
like to use two wave functions u n ( o) and such that the corresponding operator O a (uj) 

couples only to the ground and first excited state in the pseudo-scalar channel, respec¬ 
tively. As one can easily see from eq. ( |2.13| ) the sensitivity to ca is then proportional to 
m^(i) — 777^( 0) . Higher excited states are (by definition) not contributing and in principle 
the method can be used for rather small T. Hence, we find this the most attractive method 
both from a theoretical and practical point of view. In the next section we will detail our 
approximation to this ideal situation. 


2.3 Wave functions 

We will now proceed to the more technical aspects of our method. In order to approximate 
07 ^( 0 ) and u p) consider a set of N wave functions. Given a vector u in this IV-dimensional 
space, projected correlation functions are defined as (u, /a) and (u, f\u ), i.e. /a is regarded 
as a vector and f\ as a matrix in this space. It is useful to represent /x (X = A, P) and f\ 
as [31] 

M—l 

/x(x 0 ; Ui) = Y + 0{e- m ^ )xo ) + O( e ~ m G(T-x 0 )) , ( 2 . 14 ) 

n=0 

M—l 

fi(uuu;,) = Y v t ) vf ) e- m * )T + 0(e~ miM)T ) + 0(e~ m ^ T ) , (2.15) 

77.—0 

(n) 

where n labels the states in the pseudo-scalar channel in increasing energy and u) ' is the 
overlap of such a state with the one generated by the action of O a (uJi) on the SF boundary 
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state. The mass viq belongs to the lowest excitation in the scalar channel, the 0 ++ glueball 

(n) 

and the coefficients ' are proportional to the decay constant of the nth state. Here we 
have suppressed the explicit volume dependence of all quantities. 

Knowledge of v^ would allow the construction of vectors u^ n \ such that - up to correc¬ 
tions of order e~ mn T - the correlation (vS- n \ /a) receives contribution from the nth state 
only. These u may be computed from the v^ by a Gram-Schmidt orthonormalization. 
Clearly, and u ^ can then be used to approximate ui n ( o) and w n (i) ■ 

An approximation to the v^ n ' 1 can be obtained from the eigenvectors of the positive 
symmetric matrix f\. For the normalized eigenvectors 77 !°), rj^\ ... corresponding to eigen¬ 
values A(°) > A^ > .. . eq. ( |2.15| ) implies that 


|£,(°) _ ^ 


(0) 112 _ Qj' e -(ml- 1) -ml 0) )r^ 


and {rP,v^) = 0{e- { ^ ) -^ T ) . 


(2.16) 

(2.17) 


Thus, to the order indicated above, is given by r/°) and 7 / 1 ) is orthogonal to the 
’’ground state vector” As eigenvectors of a symmetric matrix the 7 /”) are already 

orthogonal and we therefore use the approximation 

w 7r(°)-XX° )wi and ( 2 - 18 ) 

i i 

to obtain correlators, which are (for intermediate xo) dominated by the ground and first 
excited state, respectively. We note in passing that the ratios r uf^ /v "' have a continuum 
limit if the wave functions are properly scaled with the lattice spacing. 

In our simulations we restrict ourselves to a basis consisting of three (spatially periodic) 
wave functions defined by 


W;(x) = N i 1 ^2 w *(l x “ nL \) > * = 1, - - -, 3 , 

ngZ 3 * 

uli(r) = 7y( 3 / 2 e” r / a ° , ^(t - ) = r^ 5 ^ 2 r e~ r ^ a ° , 

aj 3 ( r ) = r “ 3 / 2 e -”/( 2 “o) , (2.19) 


where ao is some physical length scale. We thus keep it fixed in units of L, choosing 
ao = L/6. The (dimensionless) coefficients N{ are fixed to normalize the wave function via 
a 3 D) x w?(x) = 1. In addition we also consider the flat wave function cuo(x) = L~~ 3 / 2 , where 
both quarks are projected to zero momentum separately . 1 


3. Numerical computation 

3.1 Results for the improvement coefficient 


All our simulations were performed using non-perturbatively improved Wilson fermions [5, 
7,13,14] and the plaquette gauge acti on. For the boundary-improvement coefficients ct and 

1 Since in this case x and y in eq. (2.11) are uncorrelated, full translational invariance can be used 

without performing additional inversions of the Dirac operator. For an. ..3 we replace one of the spatial 

sums in eq. (2.11) by a sum over eight far separated points, which means that one performs eight times as 

many inversions. In a dynamical fermion computation this additional effort is still small compared to the 

effort invested into the “updating”. 
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c.t we used the 2-loop [32] and 1-loop [33] values, respectively. Concerning the algorithm, 
we employed the Hybrid Monte Carlo with two pseudo-fermion fields as proposed in [34], 
For all observables we have checked the expected scaling of the statistical error with the 
sample size and thus verified the absence of the problems described in [35] at the volumes 
and masses we consider here. Table |l] summarizes the parameters of our simulations. 


run 

L/a 

T/a 

P 

K 

-^meas 

am/t({3 ) 

-c A 

i 

12 

12 

5.20 

0.135600 

320 

0.0151(9) 

0.0638(23) 

ii 

16 

16 

5.42 

0.136300 

200 

0.0171(5) 

0.0420(21) 

hi 

24 

24 

5.70 

0.136490 

120 

0.0151(4) 

0.0243(36) 

IV 

12 

12 

5.20 

0.135050 

160 

0.0363(6) 

0.0697(31) 

V 

16 

20 

5.57 

0.136496 

290 

0.0154(4) 

0.0366(36) 

VI 

24 

24 

6.12 

0.136139 

40 

0.0002(4) 

0.0244(21) 


Table 1: Summary of simulation parameters and results for ca- Runs I-III are at constant physics. 


The (3 values for run II and III have been chosen such that L/r o is approximately the 
same as in run I, which corresponds to L ~ 1.2 fm. In exploratory quenched studies [12] 
this volume was found to be sufficient for the described projection method to work. -/V meas 
is the number of estimates of ca, separated by 4-12 unit length HMC trajecrories. The 
autocorrelation of these measurements turned out to be negligible. The column labeled 


am/t(/3 ) refers to the bare quark mass m = r(T/ 2;cjo) + ocas(T/2;cjq), cf. eqs. (2.12 


2.13| ). The 1-loop value of ca from [6] is used there. We tuned the hopping parameter k in 
order to keep am/t((3 ) fixed when varying /3, thus ignoring (presumably small) changes in 
the renormalization factors. Note that we have chosen a finite, but small bare quark mass 
of around 30 MeV. Such a mass helps (in addition to the Dirichlet boundary conditions) 
to reduce the cost of the simulations. Results from the remaining simulations are used to 
discuss systematic uncertainties in our determination of c\. 

In Fig. |2]we show the effective masses from fp(xo]UJ n (o)) and /p(xo; w^p)) as obtained 
in run II. Two distinct signals are clearly visible, which indicates that the described approx¬ 
imate projection method works well at these parameters. The energy of the first excited 
state is not far away from a -1 , suggesting that in even smaller volumes the residual 0(a 2 ) 


effects would grow rapidly. In the spirit of the remark after eq. (2.18) at the other values of 
/3 we used the same linear combination of wave functions to define u; ( 0 ) and u p), namely 


rf °) = (0.5172, 0.6023, 0.6081) 

and rp) = (0.8545, -0.3233, -0.4066) 


(3.1) 


which are the ones determined in run II. When scaled in units of ro, this yields effective 
masses similar to those shown in Fig. Results from a redetermination of and t/ 1 ) 
in the other matched simulations would differ from eq. e only by ~ 1%. In the figure 
the error on the effective mass of the first excited state is seen to be quite large, but what 
actually enters the computation of ca is the error of 


Ar(x 0 ) = r(x 0 ; ) - r (x 0 ; lu ^ 0) ) and A s(x 0 ) = s(x 0 ; u nW ) - s(x 0 ; w l(0) ) . (3.2) 
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xo/a 

Figure 2: The effective mass in lattice 
units of the projected correlation functions 
(h ( °\/p) and (? 7 (1) ,/p) from run II. 



Figure 3: As(xo) and Ar(xo) determined 
from L 0 n ( o) and ui w (i) in run II. 


These profit from statistical correlations of the correlation functions entering their definition 
and thus have smaller statistical errors as can be seen in Fig. [|, where we plot aAr and 
a 2 As from the same data used in Fig. |2|. 

Fig. |4| collects results for the ’’effective” ca(£o) = — Ar(xo)/aAs(a:o) from the matched 
runs I-III. We see little variation for xq > 6a, which we take as another signal that high 
energy states which could contribute large 0(a) ambiguities in the improvement condition 
are reasonably suppressed in this region. We complete our definition of c\ with the choice 
xq = T/2, which is at the same time scaled in physical units and in agreement with the 
xo > 6 a bound for all our lattices. 

Finally, ca is plotted as a function of in Fig. [|]. The solid line is a smooth in¬ 
terpolation of the data from the matched simulations, constrained in addition by 1-loop 
perturbation theory: 

c A «) = —0.00756 (3.3) 

It is our final result, valid in the range 0.98 < < 1.16 within the errors of the data points 

(at most 0.004). 

The non-perturbative result is quite far away form 1-loop perturbation theory, which 
takes a value of ca = —0.0087 at g$ = 1.15. Using data from [36], we see that the effect on 
the result for the pseudo-scalar decay constant at this lattice spacing is as large as 10%. 

3.2 Uncertainties due to deviations from the “constant physics” condition 

We should check whether the volumes in our runs I-III are scaled sufficiently precisely or if 
systematic errors need to be added to the statistical ones on ca to cover possible violations 
of the constant physics condition. Table 1 shows that the bare PCAC mass has been kept 
constant to within about 10%. A renormalized quark mass differs by the multiplication 
with a Z-factor which is a slowly varying function of a. As explained in the beginning of 






























sect. 2, such a factor is irrelevant. Run IV is done with a quark mass which is more than 
a factor 2 larger than the one in run I, with otherwise identical parameters. The small 
difference in c\ confirms that the small deviations from the “constant mass” condition can 
be neglected. 

The other issue is the uncertainty due to 
our perturbative (or asymptotic) scaling of the 
physical length scales. It has been argued in 
sect. 2.1 that the difference to a proper non- 
perturbative scaling is rather small. Also, esti¬ 
mating a possible change by comparing 3-loop 
to 2-loop and non-perturbative scaling, gives 
a deviation in t(j3) which is smaller than 10% 
in the whole range of Fig. |l] and thus the same 
maximum deviation applies to L/a. Again this 
can be neglected altogether. Additional confir¬ 
mation comes from a comparison of the result 
from run V to our fit formula. In run V, L/a is 
20% lower than the proper value, but c\ does 
not differ significantly from the fit curve. 

Finally, by run VI we verify that the de¬ 
pendence of ca on the kinematic parameters 
disappears quickly when going to even larger values of (3. In this run we used gauge 
configurations from the calculation of Z-p [37]. Although those were produced at m = 0, 
6 = 0.5 and a much smaller volume, the resulting ca is only approximately two standard 
deviations away from our fit. 



Figure 4: Effective ca as determined from 
ui v ( o) and for runs I-III. Points with 

xo/a = 6 are marked by filled symbols. 


- 0.02 
- 0.03 
- 0.04 
$ - 0.05 
- 0.06 
- 0.07 
- 0.08 


Figure 5: Simulation results for ca- The solid line represents a fit of the data points at constant 
physics (filled circles). 
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4. Discussion 


For the 0(a)-improved action with non-perturbative c sw [13], we have determined the 
improvement coefficient c\ for (3 > 5.2, which roughly corresponds to a < 0.1 fm. The 
improvement condition was implemented at constant physics, which is necessary in the 
situation when 0(a) ambiguities in the improvement coefficients are not negligible. This 
is indeed the case here: at fi = 5.42 we have evaluated ca also for an L/a = 12, T/a = 16 
geometry. The value of ca is about a (statistically significant) 30% larger in magnitude 
than for the constant physics condition L/a = 16 = T/a. Thus, imposing the constant 
physics condition is important in the present case. 

In addition, in order to safely exclude large 0(a) ambiguities, improvement conditions 
should only involve states with energy E <C a^ 1 . On this requirement we had to compro¬ 
mise more than we would have liked to do. Our maximum values for Ea are about 0.7. 
Although this could have been improved by going to somewhat larger values of L (and T), 
this would have made the numerical computation much more expensive. 

We finally note that large 0(a 2 ) effects have been found in the Nf = 2, 0(a)-improved 
theory [17] at (3 = 5.2. These may well be related to the not so small O(a) ambiguity in 
ca that we just mentioned. This can only be investigated further by studying the scaling 
violations in quantities such as after improvement. Of course also the renormalization 
constant Za has to be known to carry out such a study. We are presently computing Za 
following the strategies of [38,39]. As an immediate application one can then 0(a)-improve 
and renormalize the bare pseudoscalar decay constants computed in [19,36,40]. 

Clearly, the method employed in this paper may also be useful to compute ca in the 
three flavor case, where c sw is known non-perturbatively with plaquette and Iwasaki gauge 
actions [14,41,42], 
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